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Abstract 

We propose a variant of the Simulated Annealing method for optimization in the multivari- 
ate analysis of differentiable functions. The method uses global actualizations via the Hybrid 
Monte Carlo algorithm in their generalized version for the proposal of new configurations. 
We show how this choice can improve upon the performance of simulated annealing methods 
(mainly when the number of variables is large) by allowing a more effective searching scheme 
and a faster annealing schedule. 

KEYWORDS: Simulated Annealing. Hybrid Monte Carlo. Multivariate minimization. 



1 



1 Introduction 



An important class of problems can be formulated as the search of the absolute minimum of a 
function of a large number of variables. These problems include applications in different fields such 
as Physics, Chemistry, Biology, Economy, Computer Design, Image processing, etc. 1 . Although 
in some occasions, such as the NP-complete class of problems |2j, it is known that no algorithm can 
surely find the absolute minimum in a polynomial time with the number of variables, some very 
successful heuristic algorithms have been developed. Amongst those, the Simulated Annealing 
(SA) method of Kirkpatrick, Gelatt and VecchijS], has proven to be very successful in a broad 
class of situations. The problem can be precisely defined as finding the value of the iV-dimensional 
vector x = (xi,X2, ...,Xn), which is an absolute minimum of the real function E(pt). For large N, 
a direct search method is not effective due to the large configuration space available. Moreover, 
more sophisticated methods, such as downhill simplex or those using the gradient of i5(x)j3], are 
likely to get stuck in local minima and, hence, might not able to reach the absolute minimum. 

SA is one of the most effective methods devised to overcome these difficulties. It allows escaping 
from local minima through tunnelling and also by accepting higher values of E(x) with a carefully 
chosen probability The method is based on an analogy with Statistical Physics: the set of 
variables (x\, . . . , xn) form the phase space of a fictitious physical system. The function E(x) is 
considered to be the system's energy and the problem is reduced to that of finding the ground state 
configuration of the system. It is known that if a system is heated to a very high temperature T 
and then it is slowly cooled down to the absolute zero (a process known as annealing), the system 
will find itself in the ground state. The cooling rate must be slow enough in order to avoid getting 
trapped in some metastable state. At temperature T, the probability of being on a state with 
energy -E(x) is given by the Gibbs factor: 

P(x) ocexp(-£(x)/T). (1) 

From this relation we can see that high energy states can appear with a finite probability at high T. 
If the temperature is lowered, those high energy states become less probable and, as T — » 0, only 
the states near the minimum of E(x) have a non-vanishing probability to appear. In this way, by 
appropriately decreasing the temperature we can arrive, when T — > 0, to the (absolute) minimum 
energy state. In practice, the method proceeds as follows: at each annealing step k there is a well 
defined temperature T{k) and the system is let to evolve long enough such that it thermalizes at 
temperature T(k). The temperature is then lowered according to a given annealing schedule T(k) 
and the process is repeated until the temperature reaches T = 0. 
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To completely specify the SA method, one should give a way of generating representative 
configurations at temperature T, and also the variation of the temperature with annealing step, 
T(k). For the generation of the configurations, the Monte Carlo method (MC) is widely used|SJ|Hl 
E| . MC introduces an stochastic dynamics in the system by proposing configuration changes x — ► x' 
with probability density function (pdf) g(x'|x), i.e, if the system variables adopt presently the value 
x, the probability that the new proposed value lies in the interval (x',x' +dx') is g(x'|x)dx'. This 
proposal is accepted with a probability /i(x'|x). Much freedom is allowed in the choice of the 
proposal and acceptance probabilities. A sufficient condition in order to guarantee that the Gibbs 
distribution is properly sampled, is the detailed balance condition: 

ff (x'|x)/i(x'|x) exp(-£(x)/T) = .g(x|x>(x|x') exp(-2?(x')/T). (2) 

Once the proposal pdf g(x'|x) has been conveniently specified, the acceptance probability /i(x'|x) is 
given as a convenient solution of the previous detailed balance equation. Usually (see next section) 
the proposal probability g(x'|x) = <?(Ax) is a symmetric function of the difference Ax = x' — x, 
(/(Ax) = g{— Ax) and a commonly used solution to the detailed balance equation is the Metropolis 
choice: 

/i(x'|x) = min (1, cxp [- (E(x') - £(x)) /T]) , (3) 
although other solutions have been also widely used in the literature. 

The various SA methods differ essentially in the choice of the proposal probability <?( Ax) and the 
annealing schedule T(k). One can reason that the cooling schedule T(k) might not be independent 
of the proposal probability g(Ax), i.e. T(k) should be chosen consistently with the selected g(Ax) 
in such a way that the configuration space is efficiently sampled. In the next section we briefly 
review the main choices used in the literature. We mention here that most of them involve only 
the change of one single variable X4 at a time, i.e. they consist generally of small local moves. N 
of these local moves constitute what is called a Monte Carlo Step (MCS). The reason for using 
only local moves is that the acceptance probability given by is very small if all the variables 
are randomly changed at once, because the change in energy E(x') — E(x) is an extensive quantity 
that scales as the number of variables N. Hence, the acceptance probability near a minimum of 
-E(x) becomes exponentially small. Since Ax is a small quantity, the cooling schedule must be 
consequently small, because a large cooling rate would not allow the variables to thermalize at the 
given temperature. It is then conceivable that the use of a global update scheme could improve 
upon the existing methods by allowing the use of larger cooling rates. 

In this paper we investigate the effect of such a global update dynamics. Specifically, we use the 
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Hybrid Monte Carlo (HMC) algorithm |Hj for the generation of the representative configurations at 
a given temperature. By studying some examples, we show that the use of this global dynamics 
allows quite generally an exponentially decreasing cooling schedule, which is the best one can 
probably reach with other methods. Another advantage of the use of the HMC is that the number 
of evaluations of the energy function -E(x) is greatly reduced. Finally, we mention that the use of 
a generalized HMC jUJ EH allows to treat efficiently minimization problems in which the range of 
variation is different for each variable. 

The rest of the paper is organized as follows: in section II we briefly review some of the existing 
SA methods; in section III we explain how to implement Hybrid Monte Carlo in an optimization 
problem; in section IV we use some standard test functions to compare our method with previous 
ones; and in section V we end with some conclusions and outlooks. 

2 Review of Simulated Annealing Methods 

Amongst the many choices proposed in the literature, we mention the following: 

-Boltzmann Simulated Annealing (BSA) 11 : Based on a functional form derived for many phys- 
ical systems belonging to the class of Gaussian-Markovian systems, at each annealing step k the 
algorithm chooses a proposal probability given by local moves governed by a Gaussian distribution: 



The Metropolis choice © is then used for the acceptance. This choice for the proposal probability 
and the use of purely local moves imply that the annealing schedule must be particularly slow: 
T(k) = To/ ln(l + Afc), for some value of the cooling rate A. 

-Fast Simulated Annealing fFSA)[T2]: States are generated with a proposal probability that has a 
Gaussian-like peak and Lorentzian long-range tails that imply occasional long jumps in configura- 
tion space. These eventual long jumps make FSA more efficient than any algorithm based on any 
bounded variance distribution (in particular, BSA). The proposal probability at annealing step k 
is a ./V-dimensional Lorentzian distribution: 



g (Ax) ~ exp 



f |Ax| 2 l 



(4) 



2T(k) 



ff(Ax) - T(fc)(|Ax| 2 + T(fc) 2 )-^. 



(5) 
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One of the most significant consequences of this choice is that it is possible to use a cooling schedule 
inversely proportional to the annealing step k, T(k) = Tq/(1 + Afc), which is exponentially faster 
than the BSA. 



-Very Fast Simulated Reannealing fVFSR)|13|: In the basic form of this method, the change Ax 
is generated using the set of random variables y = (yi, . . . , yjy) 

= {Bi - A t )yi, (6) 

(Ai and Bi are the minimum and maximum value of the i-th dimension range). The proposal 
probability is defined as 

N 1 
5(y) = J3 2(|j/ i |+T i (fc))ln(l + l/r i (fc))- (?) 

Notice that different temperatures Tj(fc) can be in principle used for the updating of different vari- 
ables Xi. For the acceptance probability, one uses the Metropolis choice © with yet another tem- 
perature To(fc). This proposal allows the following annealing schedule: Ti(k) = 2^(0) exp(— Ajfc^), 
i = 0,1,..., N, which is not very efficient for large number of variables N. A more detailed 
description of the VFSR algorithm can be found in |13| . 

-Downhill Simplex with Annealing (DSA):4-: This method combines the Downhill Simplex (DS) 
method (which is basically a searcher for local minima) with a Metropolis like procedure for the 
acceptance. The DS samples the configuration space by proposing moves of the "simplex" . A 
simplex being a geometrical figure with N + 1 vertices in the A-dimensional phase space. The 
moves are usually reflections, expansions, and contractions. The acceptance part is implemented by 
adding logarithmically distributed random variables proportional to the temperature to the energy 
before the move and subtracting a similar random variable after the move. The move is accepted 
if the energy difference is negative. According to reference 4_ different annealing schedules T(k) 
should be used for different problems. In the implementation we have made of this method (see 
section IV) an exponential decay has been used. 



3 Hybrid Simulated Annealing 



The alternative method we propose -Hybrid Simulated Annealing (HSA)- uses the Hybrid Monte 
Carlo (HMC)jHj in their generalized version E3j to generate the representative configurations. 
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We first review the HMC method. 

In its simplest and original form, HMC introduces a set of auxiliary momenta variables p = 
(pi, . . . ,pn) and the related Hamiltonian function 7i(x, p): 

1 N 

W(x, p) = E( X1 , . . . , a:*) + - £ = £(x) + p 2 /2. (8) 

From the Gibbs factor: 

P(x,p) oc exp[-«(x,p)/T] = exp[-£(x)/T] exp[-p 2 /2T], (9) 

we deduce that, from the statistical point of view, the momenta p are nothing but a set of inde- 
pendent, Gaussian distributed, random variables of zero mean and variance equal to the system 
temperature T. There is no simple closed form for the proposal probability g(x'\x), and the pro- 
posal change x — > x' is done in the following way: first, a set of initial values for the momenta p 
are generated by using the Gaussian distribution exp[— p 2 /2T] as suggested by the equation ©; 
next, Hamilton's equations of motion, Xi = pi, pi = Fi, where -Fj(x) = —dE(x)/dxi is the "force" 
acting on the variable X4, are integrated numerically using the leap-frog algorithm with a time 
step St: 

St 2 

x[ = Xi+Stpi + — Fj(x) (10) 

Pi = JX + j[fi(x)+F<(x')] > i = l,...,N. 

The proposal x' is obtained after n iterations of the previous basic integration step. In other words: 
by numerical integration of Hamilton's equations during a "time" nSt. The value x' must now be 
accepted with a probability given by: 

/i(x'|x) =min(l,exp[-(H(x',p / )-W(x,p))/T]). (11) 

Notice that this acceptance probability uses the total Hamiltonian function 7Y(x, p) instead of 
simply the function E(x) as in the methods of last section (compare and Q) 14 . Although 
Hamilton's equations exactly conserve the energy Iri — 0, the difference AH = W(x', p') — "H(x, p) 
is not equal to zero due to the finite time step discretization errors and one has quite generally 
AH. = 0(NSt l ) for some value of I. In this way, although the mapping is a global one, i.e. all the 
variables are updated at once, it is still possible to have an acceptance probability of order unity 
by properly choosing the time step St and one can have large changes in phase space at a small cost 
in the Hamiltonian. Notice that the Hamiltonian difference AH being small, does not necessarily 
imply that AE is small and once can in principle accept moves which imply a large change in the 
energy E(x). 
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In order to generate configurations at temperature T, one still must satisfy the detailed bal- 
ance condition, equation (J2J). One can prove that sufficient requirements for this detailed balance 
condition to hold are that the mapping given by eas. Hll)|) satisfies time reversibility and area 
preserving 16 . These two properties are exactly satisfied by Hamilton's equations and are also 
kept by the leap-frog integration scheme. Under those conditions, the Gibbs distribution for 
the original variables x is properly sampled. It is possible to further generalize the HMC method 
by using more general mappings satisfying the conditions of time reversibility and area preserving. 
In reference it was shown that those conditions were satisfied by the mapping induced by n 
iterations of the following basic step: 

N St 2 N 

x\ = n + 8t^A ijPj + —Y1 AikAjkFjW, (12) 
j=i j,k=i 

St N 

Pi = ft + YE^^iW+^Wl' i = l,...,N, 

i=i 

where Ay is an arbitrary matrix. This mapping can be thought as the leap-frog numerical inte- 
gration of the following equations of motion: 

±i = ^2Ai jP j, (13) 

3 

Pi = J2 A r' F J- 

3 

An straightforward calculation shows that these equations, although not being Hamiltonian, still 
conserve energy, 7i = 0, and the main features mentioned above of the standard HMC method 
are still maintained. Convenient choices for matrix Aij are: diagonal in Fourier space (Fourier 
acceleration), or a diagonal matrix: Ay = Ajtfy. This last choice allows an effective integration 
time step Sti = AiSt different for each variable (compare with (|10|0 : 

St 2 

x'i = Xi+StiPi + ^Fifr), (14) 
Pi = Pi + ^[Ji(x)+fi(x')], i=l,...,N 

The possibility of using different time steps for each variable accounts for the fact that the range 
of variation might differ for each variable. This is the case, for instance, of Corana's function (see 
next section). 

Summing up, the HMC proceeds by generating representative configurations by using a proposal 
obtained by some of the mappings given above. This proposal must now be accepted with a 
probability given by In this paper, we have used mainly the basic mapping given by (|10|l 
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except in one case (Corana's function) in which the mapping (|14|l has been used instead. The 
temperature must then be decreased towards zero as in other SA methods. Notice that in the case 
T = the random component of the evolution (the momenta variables) in Eq. (|10|l is zero and then 
the proposal coincides with that of gradient methods. 

The HMC has been extensively used in problems of Statistical Physics [17.. For our purpose 
here, we have found that the use of the previous Hamiltonian based global update of the statistical 
system associated with the energy E(~x), allows a much more effective annealing schedule and 
searching scheme than, for instance, the Boltzmann, Fast annealing and Very Fast Reannealing 
methods mentioned above. In particular we have been able to use quite generally an exponential 
annealing schedule: T(k) = Toe . Moreover, since in HMC the acceptance decision is taken 
after all the N variables have been updated, the number of energy function evaluations is greatly 
reduced. This turns out to be important in those problems in which the calculation of the energy 
function E{x) takes comparatively a large amount of computer time. 



In order to compare our algorithm with the different ones proposed in the literature, we have used 
a set of five test functions: a multidimensional paraboloid, a function from De Jong's test |18|. 
Corana's highly multi-modal functional an d two other functions with many local minima. We 
now define and describe in some detail these functions. 

The first function, /i(x), is a ^-dimensional paraboloid: 



Here we use the test value N — 200 and to compare with the results in j2D], we also use the value 
N = 3. Although this is a particularly simple function with a single minimum f± = located at 
Xi = 0, i = 1 . . . N, it ultimately describes the late stages of the behaviour of the SA algorithm 
when we are near a local or global minimum of any differentiable function. 

The second function, /2(x), is a two dimensional (N = 2) function taken from De Jong's test 



4 Results 




(15) 



typically used for benchmarking Genetic Algorithms ^B] : 



/ 2 (x) 



25 



0.002 + Y,U + («i ~ a jf + fca - bj) 6 }- 1 

3=1 



(16) 



where the vectors a, b have the following 25 components: 

a 3 = {-32,-16,0,16,32,-32,-16,0, 16, 32,..., -32, -16, 0,16, 32}, 

bj = {-32, -32, -32, -32, -32, -16, -16, -16, -16, -16, . . . , 32, 32, 32, 32, 32}, 

this function has 25 local minima, and the global minimum is fi = 0.998004, at x\ = X2 = —32. 



The /a(x) function is the Corana's function: 



N 



/3(x) = I ri ./^ 

l — l I ^iJj, 



0.15 (0.05 sgn(zi) + Zi ) 2 dk if \x> - z-\ < 0.05 



(17) 



otherwise 

z t = 0.2 LI 5 Xi\ + 0.49999jsgn(xi) (18) 

di is an iV-dimensional vector. In our tests (and following [23| we have used N = 10 and 
d = (1,1000,10,100,1,1000,10,100,1,1000). This function, which has many local minima and 
is discontinuous and piecewise differentiable, turns out to be one of the most difficult test func- 
tions, because the different variables have different scales of variation. The global minimum is 
/3(x)=0,atasi = 0,i = l...iV. 



The /4(x) function is defined by: 



/4(x) = ^y__! 



1 x—^ sin(47ri4rxj) 



(19) 



2N f-f sin(27TXi) 

with TV = 200, K = 2. This function is periodic and has (2K - l) N local minima per period. 
The absolute minima are at Xi = (2m + l)/2, m £ Z, i = 1...N, and the minimum value is 
/4(x) = — K (sec figure 0. 



And, finally, the /s(x) function is defined by: 



N 



N 



/ B (x) = J2 w Q - ll cos ( 47r ^)' ( 2 °) 

»=1 i=l 

with iV = 10 and a = 1.3. Again, this function has many local minima. The absolute minimum is 
/ 5 = -1 at Xi = 0, i = 1...N. 
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We present results of the optimization of these typical test functions performed with the meth- 
ods described above: Fast Simulated Annealing (FSA) , Very Fast Simulated Reannealing (VFSR) , 
Downhill Simplex with annealing (DSA) and the Hybrid Simulated Annealing (HSA). Amongst 
other quantities, we have focused, as usual in this field, on the number of evaluations of the function 
and the CPU time needed to achieve a given accuracy e in the minimum value of each function. 
These minimum values being exactly known for the test functions used. The results are summa- 
rized in tables and J5J) after averaging over 10 realizations. An accuracy value of e = 10~ 3 has 
been used, although similar results hold for other values of e. We have programmed the algorithms 
for the FSA, DSA and HSA methods, whereas the results for VFSR have been taken directly from 
|20|. For a given test function, we have used the same initial condition, Timitiai, for each method. 
As a general trend, we can see that HSA performs better than the other methods when the number 
of variables N is large. This does not imply that HSA performs extremely worse for small values 
of N. An important advantage of HSA in front of other methods is that the number of function 
evaluations is much smaller (in table the number of function evaluations includes also the cal- 
culation of the forces necessary in the HSA method) . This might turn out to be very important 
in those problems in which the function evaluation takes a long computer time. We now report in 
some detail the results of each test function: 

As mentioned before, the fa function, a parabolic function with a single minimum, serves 
to model the behavior close to a minimum of any function, i.e. the situation for low enough 
temperature. When the number of variables is small, N = 3, it turns out that the fastest method 
(in the sense that it reaches the minimum in less computer time) is DSA although HSA needs less 
function evaluations. However, when the number of variables is large, N = 200, the cost in CPU 
time and number of function evaluations is very favorable to HSA. In general, the performance of 
the DSA method worsens when the problem has many minima. This is obvious when looking at 
the results for the De Jong's fa, the Corana fa and the fa functions for which the DSA could not 
even find the absolute minimum. 

The fa function is another example in which the HSA can not offer a better alternative than 
other methods, stressing the fact again that for small number of variables the use of a global actu- 
alization turns out to be irrelevant. In this case, VFSR needs less number of function evaluations 
than any other method. However, for large number of variables N, the cooling schedule required 
for VFSR is necessarily slow (see the discussion in section 2) making it inefficient for large N. 

The functions in which the variables have a wide range of variation (for instance Corana's 
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function fz) can be better handled using the general version of HSA. Remember that the rescaling 
in 1)14(1 allows an effective integration time step SU = AiSt different for each variable. So, one can 
tune Aj to solve efficiently this kind of problems. In our case, the range of variation of the variables 
come essentially from the part V(x) = dixf of the Corana's function. Then, from the equations 
of motion £4 = Aipi and p\ — AiFi we have Xi = AfFi. The force is Fi — —2 diXi and we have 
ii — —AfdiXi , so we chose Ai — -h=- in order that each variable has the same effective time scale 
for evolution. 

The fi and f$ functions have the feature of possessing a large number of minima (for example, 
/4 has (2K — 1) N local minima in a period). The results show again that HSA is a much better 
alternative when the number of variables is large, both from the point of view of CPU time used or 
the number of function evaluations. We have chosen the function to compare in figure J2J the 
evolution of the minimum value of the function with the actual number of function evaluations, 
for both the FSA and HSA methods, showing again in a different manner that HSA can find a 
better minimum with a less number of function evaluations. From the results for these functions 
we infer that in minimization problems with a large number of variables and a large number of 
local minima, the HSA has the best performance. Needless to say, we have made our best effort 
to use the optimal values for the parameters in each method. It is possible, though, that these 
values could be further improved and the results of tables (JIJ and J5J slightly modified. We believe, 
though, that this will not affect the main conclusions of this paper. 



5 Conclusions 



We have shown by some examples how the use of the global update using Hybrid Monte Carlo 
algorithm can indeed improve the performance of simulated annealing methods. The global updat- 
ing implicit in HSA allows an effective searching scheme and fast annealing schedules and becomes 
highly effective, mainly in those problems with a large number of variables and a large number of 
metastable minima. 

It is clear from the results in the previous section that HSA requires in some cases orders of 
magnitude less evaluations of the function than other methods and can, therefore, give a solutions 
in less computer time. This conclusion remains despite the fact that HSA requires some extra work 
when computing the evolution equations since it needs to compute also the forces Fi acting on 
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the different variables. In those cases in which the evaluation of the function takes a considerable 
amount of computer time, HSA will have an optimal performance, since the number of function 
evaluations is greatly reduced as compared to other simulated annealing methods. It is conceivable 
also that one could then use efficiently some of the acceleration schemes (Fourier, wavelet, etc.) 
available for Monte Carlo methods in order to improve upon the convergence of the simulated 
annealing techniques. Further developments include applying HSA to techniques such as the 
Car-Parrinello method for finding the ground state of quantum many body systems, for which 
the calculation of the energy function is very time consuming. Work on this direction is under 
progress. 
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Function 


Dimension 


FSA 


VFSR 


DSA 


HSA 


fx 


3 


480 


4875 


79 


18 


fx 


200 


8420000 




474000 


30 


h 


2 


9900 


1476 


(*) 


165000 


h 


10 


2100000 


319483 


(*) 


720000 


h 


200 


12925000 




(*) 


163000 


h 


10 


7230000 




570000 


118000 



Table 1: Number of function evaluations averaged over 10 realizations, for each of the simulated 
annealing methods used for optimization of the different functions to reach the absolute minimum 
with an accuracy of e = 10~ 3 . In those cases marked (*) it was not possible to reach the absolute 
minimum. For the HSA, the displayed number is the number of function evaluations including the 
calculations of the force. 



Function 


Dimension 


FSA 


DSA 


HSA 


fx 


3 


0.023 


0.003 


0.021 


fx 


200 


182.898 


163.763 


0.039 


h 


2 


0.181 


(*) 


5.834 


h 


10 


29.454 


(*) 


11.730 


h 


200 


1662.177 


(*) 


61.863 


h 


10 


119.434 


13.929 


4.00 



Table 2: Similar to table Q but showing the CPU time (in seconds) needed to reach the absolute 
minimum with an accuracy e = 10~ 3 for each of the simulated annealing methods explained in the 
text. All the programs were run on a Silicon Graphics Origin200 (CPU: R10000 running at 180 
MHz, Speed: 15.5 SPECfp95). 
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Figure 1: Plot of fa function, equation (|19fl for N — 2, in one period. Notice the presence of many 
relative minima, but only one absolute minimum at x\ = xi = 0.5. 
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1 0.000 E 




Number of Function Evaluations 

Figure 2: Plot of "energy" difference with respect to the ground state value, versus the number 
of function evaluations, for the function l|19|l with N — 200 using HSA (dotted line) and FSA 
(continuous line), both initialized in Xj„jtj a ; = 1.0, the other parameters have the following values, 
for FSA: T = 0.8, m = 100, A = 100; and for HSA: T = 1.0, m = 10, n = 10, St = 0.3 A = 0.007, 
where m is the number of MCS used for thermalization at temperature T(k) 
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